ChIP AML PiPeline v2¶
import os
import pandas as pd
import sys
sys.path.insert(0, '../..')
import itertools
from scipy import stats
import numpy as np
from genepy.epigenetics import chipseq as chip
from genepy.utils import helper
from genepy.utils import plot
from genepy.epigenetics import plot as cplot
import igv
import SimpSOM as sps
from scipy import stats
import seaborn as sns
from matplotlib import cm
from matplotlib import pyplot as plt
from IPython.display import IFrame
import seaborn as sns
from bokeh.plotting import *
import igv
import numba
from numba import jit
from scipy.cluster.hierarchy import linkage, leaves_list
from sklearn.cluster import AgglomerativeClustering, DBSCAN, KMeans, OPTICS
from sklearn.mixture import GaussianMixture
from sklearn.manifold import MDS, TSNE
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
from IPython.display import IFrame
from statsmodels.stats.multitest import multipletests
from pybedtools import BedTool
import pyBigWig
output_notebook()
%load_ext autoreload
%autoreload 2
#setting basic parameters
project="Cobinding_ChIP"
version="v3"
merging_version="remove_single_wmotif"
window="150"
merged = pd.read_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/merged_with_motifs.bed.gz', sep='\t', index_col=0)
merged.columns
#retrieving the data from the previous notebook
%store -r merged
%store -r chrombed
%store -r mergedpeak
%store -r cols
%store -r annot
%store -r version
%store -r merging_version
%store -r window
%store -r crc
merging_version
merged = merged.fillna(0)
# get onmotif regions
rename = {}
torm= []
for v in merged.columns:
if "_wmot" in v:
merged[v.replace('_wmot','')] = merged[v]
rename.update({v.replace('_wmot',''):v})
torm.append(v.replace('_wmot','_mot'))
merged = merged.drop(columns=rename.values()).rename(columns=rename).drop(columns=torm)
merged[merged.columns[cols:]] = merged[merged.columns[cols:]].astype(float)
merged = merged[merged[merged.columns[cols:]].astype(bool).sum(1)!=0]
#merged = merged[merged[merged.columns[cols:annot]].astype(bool).sum(1)>1]
merged = merged.reset_index(drop=True)
helper.createFoldersFor("../results/"+project+"/"+version+"_"+merging_version+"_"+window+"/plots/")
Correlation with Annotations¶
non binarized¶
#data = pd.DataFrame(np.corrcoef(stats.zscore(0.01+merged[merged.columns[cols:]].T)), columns=merged.columns[cols:], index = merged.columns[cols:])
data = pd.DataFrame(np.corrcoef(merged[merged.columns[cols:]].T.values.astype(float)), columns=merged.columns[cols:], index = merged.columns[cols:])
link = linkage(data[:annot-cols])
col = data.iloc[annot-cols:]
col = col[[co for co in col.columns if co not in col.index.tolist()]]
rdb = cm.get_cmap('RdBu_r', 256)
for val in col.columns:
a = [rdb(128+int(v*128)) for v in col[val]]
col[val] = a
Correlation: hierarchical clustering¶
A hierarchically clustered heatmap representing the correlation coefficient of each CRC TF's signal with each other within the cobinding matrix. Above is displayed the correlation of non CRC proteins/marks with CRC TFs.
#correlation_withannotation
fig = sns.clustermap(data.iloc[:annot][data.columns[leaves_list(link)]], vmin=-1,vmax=1, col_colors=col.T, figsize=(10,20), cmap='RdBu_r', row_linkage=link, col_cluster=False, colors_ratio=0.01)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("correlation_withannotation")
#fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/correlation_with_annotation.pdf')
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/boolcorrelation_with_annotation.pdf')
plt.show()
binarized¶
similar plot computed on binarized signal (whether a peak is present at a location or not)
#data = pd.DataFrame(np.corrcoef(stats.zscore(0.01+merged[merged.columns[cols:]].T)), columns=merged.columns[cols:], index = merged.columns[cols:])
data = pd.DataFrame(np.corrcoef(merged[merged.columns[cols:]].T.astype(bool)), columns=merged.columns[cols:], index = merged.columns[cols:])
link = linkage(data[:annot-cols])
col = data.iloc[annot-cols:]
col = col[[co for co in col.columns if co not in col.index.tolist()]]
rdb = cm.get_cmap('RdBu_r', 256)
for val in col.columns:
a = [rdb(128+int(v*128)) for v in col[val]]
col[val] = a
#correlation_withannotation
fig = sns.clustermap(data.iloc[:annot][data.columns[leaves_list(link)]], vmin=-1,vmax=1, col_colors=col.T, figsize=(10,20), cmap='RdBu_r', row_linkage=link, col_cluster=False, colors_ratio=0.01)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("correlation_onbinarized_signal_withannotation")
#fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/correlation_with_annotation.pdf')
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/boolcorrelation_onbinarized_signal_with_annotation.pdf')
plt.show()
Looking at peak overlap¶
How many of peak in A (column) overlaps with peak in B (rows)
overlap, correlation = await chip.pairwiseOverlap(merged, norm=True)
enrichment, pvals = await chip.enrichment(merged)
#saving the enrichments
overlap.to_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/pairwise_overlap.csv')
correlation.to_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/pairwise_correlation_on_overlap.csv')
enrichment.to_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/pairwise_enrichment_on_overlap.csv')
pvals.to_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/pairwise_enrichment_on_overlap_pvals.csv')
merging_version = "remove_single"
overlap = pd.read_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/pairwise_overlap.csv', index_col=0)
correlation = pd.read_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/pairwise_correlation_on_overlap.csv', index_col=0)
enrichment = pd.read_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/pairwise_enrichment_on_overlap.csv', index_col=0)
pvals = pd.read_csv('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/pairwise_enrichment_on_overlap_pvals.csv', index_col=0)
pairwise overlap¶
Overlap of ROW proteins peaks in COLUMN protein peaks. A hierarchically clustered heatmap where 1 represent 100% overlap and 0, 0%. (e.g. E2F3 binds almost always (>80% of the time) under SREBF1. But SREBF1 only does so ~10% of the time.) The bottom-left side square are the CRC (clustered). At the top and right side has been added additional signals of histone marks/state and non CRC proteins (non clustered). The plot is non symmetrical as a protein can bind in many more places than another and thus cover all its binding regions while binding in many other places.
link = linkage(overlap.iloc[:annot-cols]) # D being the measurement
c = overlap[[co for co in overlap.columns if co not in overlap.iloc[:annot-cols].index.tolist()]]
viridis = cm.get_cmap('viridis', 256)
for val in c.columns:
a = [viridis(int(v*255)) for v in c[val]]
c[val] = a
fig = sns.clustermap(overlap.iloc[:annot-cols][overlap.columns[np.concatenate((leaves_list(link), range(annot, len(overlap.columns))))]], row_linkage=link, col_colors=c,figsize=(20,18), colors_ratio=0.01, col_cluster=False)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("pairwise_overlap_clustermap")
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/pairwise_overlap_clustermap.pdf')
plt.show()
Correlation on overlaps¶
Here we display the correlation on binding signal, only on regions where the two proteins/signals overlap. It allows us to recover signal that would be hidden for sets of TF that would only co-bind in few places.
link = linkage(correlation.iloc[:annot-cols], optimal_ordering=True) # D being the measurement
c = correlation[[co for co in correlation.columns if co not in correlation.iloc[:annot-cols].index.tolist()]]
viridis = cm.get_cmap('viridis', 256)
for val in c.columns:
a = [rdb(128+int(v*128)) for v in c[val]]
c[val] = a
fig = sns.clustermap(correlation.iloc[:annot-cols][correlation.columns[np.concatenate((leaves_list(link), range(annot, len(correlation.columns))))]], row_linkage=link, col_colors=c, colors_ratio=0.01, figsize=(20,18), vmin=-1, vmax=1, cmap='RdBu_r', col_cluster=False)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("correlation_on_pairwise_overlap_clustermap")
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/correlation_onoverlap.pdf')
plt.show()
Enrichments on the cobinding matrix¶
Log-enrichments given for each TF and chromatin mark over each other. Enrichment is computed by doing a fisher's exact test on the expected vs observed presence of ROW protein peaks under COLUMN protein peaks. Pvalues are then corrected for multiple hypothesis testing using the BH method. (e.g. here CTCF is very strongly enriched in SMC1 but SMC1 is not as much)
crc = crc & set(enrichment.columns)
sns.boxplot(data=enrichment.loc[crc][["super_enhancer","promoters","enhancer"]])
plt.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/average_enrichment_ctf_in_enhancers.pdf')
link = linkage(enrichment, optimal_ordering=True)
fig = sns.clustermap(enrichment,figsize=(20,20), vmin=-5,vmax=5, cmap='RdBu_r', col_linkage=link, row_linkage=link)
fig.ax_col_dendrogram.set_visible(False)
fig.ax_row_dendrogram.set_visible(False)
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/enrichment_clustermap_all_to_all.pdf')
plt.show()
interactive plot¶
(see above for information) this plot shows the same enrichments as above, except that it shows a decreasing square size for any pvalue above $10^{-9}$. hovering shows enrichment (value) and pvalue (pval) for each protein/mark pair
link = linkage(enrichment, optimal_ordering=True)
plot.correlationMatrix(enrichment.loc[enrichment.columns[leaves_list(link)]][enrichment.columns[leaves_list(link)]].values, dataIsCorr=True, names=enrichment.columns[leaves_list(link)].tolist(), pvals=pvals[enrichment.columns[leaves_list(link)]].loc[enrichment.columns[leaves_list(link)]].values, size=12, title= 'enrichment', folder='../results/'+project+'/'+version + '_' + merging_version + '_' + window+"/plots/", interactive=True, minval=-4, maxval=4)
looking only at crcs¶
enrichment = enrichment[merged.columns[cols:annot]].loc[merged.columns[cols:annot]]
link = linkage(enrichment, optimal_ordering=True)
plot.correlationMatrix(enrichment.loc[enrichment.columns[leaves_list(link)]][enrichment.columns[leaves_list(link)]].values, dataIsCorr=True, names=enrichment.columns[leaves_list(link)].tolist(), pvals=pvals[enrichment.columns[leaves_list(link)]].loc[enrichment.columns[leaves_list(link)]].values, size=12, title= 'enrichment', folder='../results/'+project+'/'+version + '_' + merging_version + '_' + window +"/plots/",interactive=True,minval=-4,maxval=4)
clustering¶
I have tried gaussian mixtures and Agglomerative clustering algorithm. Only the second can create a hierarchical clustering.
It seems that gaussian mixture makes more sense given the data we have, for now, is more "homogeneous".
I am still not so happy with the clustering. It can be because of the how much importance, outlier values and the high number of noisy values from locations with no peaks.
We can use similar methods to RNAseq to improve this (clamping values, log transform, first round of PCA..)
using KMeans¶
data= merged[merged.columns[cols:annot]].values
# using score
#scaled_data = stats.zscore(data)
# using regular scaling
scaled_data = (data-data.min(0))/(data.max(0)-data.min(0))
#binary scaling
#scaled_data = data.astype(bool)
#https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html#sklearn.cluster.KMeans
n = 3
n_clust=[10,20,50,100,300]
kmean = KMeans(n_clusters=n_clust[n], n_jobs=8)
groups = kmean.fit_predict(scaled_data)
centroid = kmean.cluster_centers_
plotting center location of clusters (kind of enrichment)¶
Heatmap of the centroid location of each cluster. centroids represent the center of mass of a cluster in kmean clustering. Here a cluster having a value of 1 means that the center of mass is completely shifted toward that particular TF. It showcase in what type of enhancers the cluster is clustering on (archetypical cobindings). This is very sensitive to the number of clusters and scaling applied
df = pd.DataFrame(centroid,columns=merged.columns[cols:annot]).T
#df = df/df.max()
fig = sns.clustermap(df, vmax=1)
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/kmeans_'+str(n_clust[n])+'_centroids.pdf')
showing the cobinding matrix with the clusters¶
the cobinding matrix, clustered using the kmeans clusters. (see cobinding matrix description)
rand = np.random.choice(merged.index,5000)
subgroups = groups[rand]
sorting = np.argsort(subgroups)
rainb = cm.get_cmap('gist_rainbow', len(set(groups)))
colors = [rainb(i) for i in subgroups[sorting]]
viridis = cm.get_cmap('viridis', 256)
data = merged[merged.columns[annot:]]
for val in data.columns:
data[val]=np.log2(1+data[val])
data = (data-data.min())/(data.max()-data.min())
data = data.iloc[rand].iloc[sorting]
for val in data.columns:
a = [viridis(int(v*256)) for v in data[val]]
data[val] = a
data["clusters"] = colors
Computing enrichment on the clusters¶
log-enrichment of each TF/chromatin marks under each cluster (see previous enrichment plot for description) enrichment are filtered for adjusted p value below $10^{-3}$
enr, pvals = await chip.enrichment(merged, groups=groups)
showing enrichment over cobinding¶
clustered cobinding matrix where signal for non CRC peaks has been replaced by enrichment of those over the clusters and where signal for CRC peaks has been averaged over the clusters
cplot.andrew(groups, merged, annot, enr=enr, pvals=pvals, cols=8, precise=False, folder = '../results/' + project + '/' + version + '_' + merging_version + '_' + window + '/plots/kmeans_', title = "sorted clustermap of cobindings clustered with Kmeans")
maxv=4
rdb = cm.get_cmap('RdBu_r', 256)
rdy = cm.get_cmap('RdGy_r', 256)
vals = merged.columns[cols:annot]
col = enr.loc[merged.columns[annot:]].copy().T
col[col<-maxv]=-maxv
col[col>maxv]=maxv
col = col/maxv
collink = linkage(enr.loc[merged.columns[cols:annot]].T, metric="cosine")
rowlink = linkage(enr.loc[merged.columns[cols:annot]], metric="cosine")
rowlinkAdd = linkage(enr.loc[merged.columns[annot:]])
for val in merged.columns[annot:].tolist():
col[val] =[rdb(128+int(v*128)) for v in col[val]]
col = col[[col.columns[i] for i in leaves_list(rowlinkAdd)]]
fig = sns.clustermap(enr.loc[merged.columns[cols:annot]], figsize=(25,25), vmax=4, vmin=-4, cmap='RdYlGn_r', col_colors=col, colors_ratio=0.01, col_linkage=collink, row_linkage=rowlink)
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/motifs_and_crc_enrichment_on_clusters_'+str(n_clust[n])+'_v1.pdf')
plt.show()
f, ax = plt.subplots(figsize=(20,10))
sns.countplot(groups.astype(str).tolist(), ax=ax, order=leaves_list(collink).astype(str).tolist())
plt.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/custer_size_'+str(n_clust[n])+'.pdf')
Using SOMs¶
Kohonen Self-Organizing Maps. Each nodes learn to recognize a pattern in the dataset, unregarding of its distribution. Each node tries to be as similar as possible to its close neighboor and dissimilar to nodes further away. This allows to reduce dimensionality of the dataset, unregarding of distances. This makes a lot of sense for our data type. we end up with a map of 400 datapoints which represent the binding code in the dataset.
#Import the library
size = 20
#Build a network 20x20 with a weights format taken from the raw_data and activate Periodic Boundary Conditions.
net = sps.somNet(size,size, scaled_data, PBC=True)
#Train the network for 10000 epochs and with initial learning rate of 0.01.
net.train(0.01, 10000)
#Save the weights to file
net.save('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/cobinding_SOMweights_'+str(size))
net = sps.somNet(0,0, scaled_data, loadFile='../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/cobinding_SOMweights_'+str(size), PBC=True)
pathplot = "../results/"+project+'/'+version+'_'+merging_version+'_'+window+"/plots/somNets/"
! mkdir $pathplot
SOMnode enrichment¶
displays a map of each node, respecting 2D neighboorhood. and its enrichment for a specific protein/mark.
#Print a map of the network nodes and colour them according to the first feature (column number 0) of the dataset
#and then according to the distance between each node and its neighbours.
for col in range(scaled_data.shape[1]):
net.nodes_graph(colnum=col, printout=True, show=True, path="../results/"+project+'/'+version+'_'+merging_version+'_'+window+"/plots/somNets/", colname=data.columns[col])
average differential enrichment¶
displays similarity in what the nodes have learnt. how different the sum of their weight is compared to their neighboors.
diffs = net.diff_graph(show=False, returns=True)
plt.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/cobinding_SOM'+str(size)+'.pdf')
interactive somPlot¶
displays the same differential map as above. but allows the user to look at the key TF/marks that are enriched in the node (might be more than one) this is very sensitive on the filter used to defined what is enriched and what is not.
plot.SOMPlot(net, size=size, colnames=merged.columns[cols:annot], folder='../results/' + project + '/' + version + '_' + merging_version + '_'+window+"/plots/onbool_", minweight=0.4)
enrichment of SOMnodes¶
Which node is mostly enriched in which specific signal? the color represent the amount of points mostly recognized that node. the more points, the stronger the enrichment.
Here we have which node recognize signal from super enhancers
prj=np.array(net.project(scaled_data[merged['super_enhancer']!=0]))
plot.bigScatter(prj,binsize=0.5,showpoint=False,precomputed=False, logscale=True, title='density plot of enhancers in TF cobinding space with TSNE', folder='../results/'+project+'/'+version+'_'+merging_version+'_'+window+"/plots/")
#Cluster the datapoints according to the Quality Threshold algorithm.
clusts = net.cluster(scaled_data, type='KMeans',numcl=n_clust[n])
sorting = []
for clust in clusts:
sorting.extend(clust)
sorting = np.array(sorting)
viridis = cm.get_cmap('gist_rainbow', len(clusts))
colors=[]
for i, clust in enumerate(clusts):
colors.extend([viridis(i)]*len(clust))
colors = np.array(colors)
viridis = cm.get_cmap('viridis', 256)
data = merged[merged.columns[annot:]]
for val in data.columns:
data[val] = np.log2(1+data[val])
data = (data-data.min())/(data.max()-data.min())
rand = np.random.choice(range(len(sorting)),5000)
rand.sort()
data = data.iloc[sorting[rand]]
for val in data.columns:
a = [viridis(int(v*256)) for v in data[val]]
data[val] = a
data["clusters"] = [tuple(i) for i in colors[rand]]
#sorted clustermap of cobindings clustered with SOM
fig = sns.clustermap(np.log2(1.01+merged[merged.columns[cols:annot]].iloc[sorting[rand]].T), vmin=0,vmax=3,figsize=(20,15),z_score=0, colors_ratio=0.01, col_cluster=False,col_colors=data, xticklabels=False)
fig.ax_col_dendrogram.set_visible(False)
fig.fig.suptitle("sorted clustermap of cobindings clustered with SOM")
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/SOM_'+str(n_clust)+'_onbool__clustermap_cobinding.pdf')
plt.show()
#let's look at the enrichment over SOM clusters
groups = []
for i, val in enumerate(clusts):
groups.extend([i]*len(val))
enr, pvals = await chip.enrichment(merged, groups = np.array(groups))
fig = sns.clustermap(enr, figsize=(12,12), vmax=4, vmin=-4, cmap='RdBu_r')
fig.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/kmeans_'+str(n_clust[n])+'_onbool_enrichment_on_SOMcluster_metasubset.pdf')
plt.show()
Plot TSNE density map¶
rand = np.random.choice(merged.index,30000)
red_data = TSNE(2,600,verbose=10,n_iter=1500).fit_transform(scaled_data[rand])
np.save('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/red_data.npy',red_data)
# Let's look at the data. the density map of the distribution of all pseudo enhancer over their tf binding profile
_, ax = plt.subplots(figsize=(10,10))
fig = sns.kdeplot(red_data[:,0], red_data[:,1], shade=True, ax=ax)
plt.savefig('../results/'+project+'/'+version+'_'+merging_version+'_'+window+'/plots/density_TSNE_peaks.pdf')
A density plot where each point represent 2Dimensional hexagonal bin and the intensity represents the accumulaation of datapoints in these bins. The datapoints represent a random distribution of ~10% of all peaks of the cobinding matrix. They dimensionality reduced using TSNE over the scaled signal of bound TF on these peaks.
These islands and clusters represent categories / types of enhancers defined by their binding code.
# now let's look at it from anothe rploting method. we can see many many islands
plot.bigScatter(red_data,binsize=0.4,showpoint=False,precomputed=False, logscale=True, title='density plot of enhancers in TF cobinding space with TSNE', folder='../results/'+project+'/'+version+'_'+merging_version+'_'+window+"/plots/")
## TODO: color by feature
## TODO: color by dependencies/pathways/expression.. after ABC model,